arXiv:1501.06868vl [nlin.AO] 27 Jan 2015 


Partial synchronization phenomena in networks 
of identical oscillators with non-linear coupling 


Celso FreitaJ^ and Elbert Macait 
Associate Laboratory for Computing and Applied Mathematics - LAC, 
Brazilian National Institute for Space Research - INPE, Brazil 


Arkady Pikovsk 
Department of Physics and Astronomy, 

University of Potsdam, Germany and 
Department of Control Theory, Nizhni Novgorod State University, 

Gagarin Av. 23, 606950, Nizhni Novgorod, Russia 
(Dated: January 28, 2015) 

Abstract 

We study a Kuramoto-like model of coupled identical phase oscillators on a network, where 
attractive and repulsive couplings are balanced dynamically due to nonlinearity in interaction. 
Under a week force, an oscillator tends to follow the phase of its neighbors, but if an oscillator is 
compelled to follow its peers by a sufficient large number of cohesive neighbors, then it actually 
starts to act in the opposite manner, i.e. in anti-phase with the majority. Analytic results yield 
that if the repulsion parameter is small enough in comparison with the degree of the maximum 
hub, then the full synchronization state is locally stable. Numerical experiments are performed to 
explore the model beyond this threshold, where the overall cohesion is lost. We report in detail 
partially synchronous dynamical regimes, like stationary phase-locking, multistability, periodic and 
chaotic states. Via statistical analysis of different network organizations like tree, scale-free, and 
random ones, and different network sizes, we found a measure allowing one to predict relative 
abundance of partially synchronous stationary states in comparison to time-dependent ones. 
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Lead paragraph: Large population of coupled oscillators synchronizes if the 
coupling is attractive and desynchronizes if it repulsive. However, the nature 
of coupling may depend on the strength of the compelling field: if the force on 
the oscillator becomes strong, it can switch its behavior from a “conformist” to 
a “contrarian”. We study such a population on a network. Here the oscillators 
connected to many others become contrarians first, so that synchrony breaks. 
We show that the regimes of appearing partial synchrony can be rather complex, 
with a large degree of multistability, and suggest a network measure which allows 
to predict relative abundance of static and dynamic regimes. 


INTRODUCTION 


Q, 


In a seminal work Uj, aiming to understand synchronization phenomena, Kuramoto pro¬ 
posed a mathematical model of non-identical, nonlinear phase-oscillators, mutually coupled 
via a common mean field. Studying this system, he identified a synchronization transition 
to an oscillating global mode when the coupling strength is larger than a critical value, 
which is proportional to the range of the distribution of the natural frequencies. Over the 
time, subsequent outcomes based on Kuramoto propositions have shown that his approach 
can be used as a framework to several natural and technological systems where an ordered 
behavior (synchronization) emerges from the interactions of a large number of dynamical 


agents [2|, |3j. Furthermore, works have shown that the Kuramoto model can be exploited as 
a building block to develop highly efficient strategies to process information 140. 

Recently, generalizations of the Kuramoto model toward interconnections between the el¬ 
ements more complex than the mean field one, have received considerable attention. Indeed, 
in many real-world problems, each dynamical agent interacts with a subset of the whole en¬ 
semble jg- 8], which can better described using a network topology. A variety of studies 
have analyzed the onset of the synchronization regime in this context. For a general class 
of linearly coupled identical oscillators, the Master Stability Function, originally proposed 
by Pecora and Carrol [^], allows one to determine an interval of coupling strength values 
that yields complete synchronization, as a function of the eigenvalues of Laplacian matrix 
of the coupling graph. For networks of oscillators with non-identical natural frequencies, 


Jadbabaie et al. 


10] were able to give similar bounds for the coupling strength of the Ku- 
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ramoto model without the assumption of infinitely many phase-oscillators. Among related 
works, Ref. jll] deals with a model whose natural frequency oscillators change with time, 
even when they are isolated. Ref. |12| explores the effects of delay in the communication 
between oscillators. Besides, Ref. [l3| builds a bridge between graph symmetry and cluster 
synchronization. 

Taking into consideration all of these previous results, one can roughly state that the 
Kuramoto transition to synchronization always happens if the coupling between oscillators 
is attractive; while, when it changes to be repulsive, this synchronization state is absent 
14 115] . However, the structure of the coupling can nontrivially depend on the level of 


synchrony itself. Sue 
in recent theoretical 


1 a dependence, called nonlinear coupling scheme, has been explored 


16dl9l] and experimental |20j, 2l| studies dealing with setup of global 


coupling. The main effect here is the partial synchrony, which establishes at moderate 
coupling strengths, where the coupling is balanced between the attractive and repulsive one. 

Here, we consider the effects of the non-linear coupling on a network: a set of identical 
oscillators, which communicate via a connected simple coupling graph, each element is 
forced by a (local) mean field, which encompasses the oscillators that are connected to 
it. The coupling function is tailored so that its influence is attractive, if the local acting 
field is small, or repulsive, otherwise. This coupling strategy implies that only nodes with 
a large enough number of connections may become repulsive. Thus, the hubs play a key 
role for the ensemble dynamics. A non-linear coupling parameter in the system tunes the 
critical quantity of connections and how cohesive this mean field must be in order to allow 
this transition. So, our approach can be considered as a dynamical generalization of the 
inhomogeneous populations of oscillators consisting of conformists and contrarians 22]. 
However, the type of oscillator’s depends on the force acting on it. 

Overall dynamics in the model can be qualitatively understood as follows: Let us assume 
initially that all the mean fields are small. Then, there are only attractive interactions 
( conformists ) in the system. So, in a first moment, they start to mutually adjust their 
phases. Above a threshold, the most connected oscillators start to feel a repulsive effect 
that drives them away from the synchronous state. In other words, if an oscillator has a 
sufficiently large number of neighbors and if it suffers enough cohesive pressure from them, 
instead of attractiveness, it becomes a contrarian, wishing to be in anti-phase with the 
force. Then, due to the repulsiveness of some nodes, other mean fields may also become 
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smaller. Finally, this tendency can shift nodes to attractiveness again. As a consequence, 
an intermediate configuration may emerge due to the balance these conflicting tendencies in 
the system. 

Depending on the non-linear coupling parameter, we report a variety of qualitative dy¬ 
namic behaviors. In general, for small values of the non-linear coupling parameter, we 
observed full synchronization and phase-locked states. When this parameter is increased, 
multistability, periodic and chaotic dynamics take place. 

The paper is organized as follows. Initially, we discuss the basic details of the model in 
section El In section IIIII the analytical result about the stability of full synchronization is 
presented. Numerical experiments in section HVl illustrate different possible regimes that the 
present model can display. Finally, in section HV Cl we perform a numerical exploration to 
address the correlation of stationary phase locking states with partial synchronization with 
the network parameters, by exploring different network topologies and sizes. 

II. MODEL OF OSCILLATOR NETWORK WITH NONLINEAR COUPLING 

Let us consider a system of N identical phase-oscillators represented by (0 1; ..., d N ) e T N 
coupled through a simple and connected undirected graph A. The dynamics for the i-th 
oscillator is given by the following ordinary differential equation 



( 1 ) 


jeMi 


where Mi denotes the set of neighbors of i in the coupling graph A. Equations (P) are 
formulated in the reference frame rotating with the common frequency of the oscillators, so 
that the latter one does not appear in the equations. The time is normalized by the linear 
coupling strength. The main feature of our model (JT]) is the nonlinear coupling parameter , 
£ > 0, which modifies the coupling at each node, with £ = 0 a standard setup of the 
Kuramoto model on a network is recovered. This parameter enters as a coefficient at the 
square of the norm of the local order parameter: 



( 2 ) 


j'eM 
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which measures the magnitude of the force acting on oscillator with index i. On the other 
hand, we represent the (global) order parameter by 


1 N 

ffe^ = 4Ve ie ' 
N ^ 


( 3 ) 


i=1 


where R G [0,1] is its norm and ip G (0, 2n] is its phase. 

We stress that Zi is not normalized (in the sense that there is no division by the number 
of the terms in the summation, like in R ), as it measures the total action of the neighbors on 
the i-th oscillator, which is called local mean field. Simple calculations show that Zf G [0, d 2 ], 
where di denotes the degree of the i-th vertex, that is, the number of incoming or outgoing 
connection, since the graph is undirected. Thus, a necessary condition for a node to suffer 
repulsive coupling, i.e. (1 — eZf) 2, < 0, is that e > df 2 . 

The introduced order parameters R and Z l} , Z n are maximal in the case of full syn¬ 
chronization 0i — ... — 9 n, whiles they decrease when oscillators begin to move apart from 
each other. 

If x < 1, where Z( nax := ma x{Zf ,..., Z'f ; }, then all oscillator will attract each other 
so that the full synchronization is established. Next, if Zf nax becomes larger than (e) _1 , the 
corresponding oscillator begins to be repulsive related to its local mean held, and the full 
synchronization breaks. As a result, may decrease and switch again the node to be 

attractive. Depending on the coupling graph A, on the initial conditions ( 9 (,..., and 
on the intensity of the nonlinear coupling parameter e, numerical simulation reveals that 
model ([TJ) can exhibit different qualitative behaviors. 

If the largest degree in the coupling graph satisfies e < d ~ with d max := maxjdi,..., djy}, 
then no node can be repulsive. We demonstrate in the next section via the Lyapunov anal¬ 
yses, that in fact this condition guarantees that the full synchronized state is stable. 


III. STABILITY OF FULL SYNCHRONIZATION 

The basic steps to obtain the results in this section follow the main ideas from [jdj. We 
begin presenting some preliminary concepts, including elements of the graph theory needed, 
and a generalized norm of the order parameter to define our Lyapunov function. 

Let B be the directed incidence matrix of a graph A. Thus, B is a matrix with N rows and 
E columns, where E is the number of directed edges of the matrix. The number of undirected 
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FIG. 1. Example of graph with N = 3 and E = 4, its directed incidence matrix B\ and its 
Laplacian matrix L\. 


edges , i.e., ignoring the direction, equals is E/2. The columns of B represent the edges of 
the graph: if the fc-th arrow (directed edge) of the graph goes from i to j, then the k -th 
column of B is zero, except at positions i and j, where B& = 1 and B ]k = —1. Regarding 
the dynamics of the system, an arrow from node i to node j in the graph means that node 
i influences node j. Although the directed incidence matrix is generally defined for directed 
graphs, it must be emphasized that only undirected graphs are considered here. We abuse 
terminology and identify a graph A with its adjacency matrix , which is an N X N matrix 
where An = 0; A %3 = A 3i = 1, if there is and edge between nodes i,j ; and A tJ = A 3i = 0, 
otherwise. So, E = J2fj =i Aj • Another common characterization of a graph is its Laplacian 
matrix , L := diag(di,..., d n) — A. One can check that L = l /zBB T . A simple illustration 
of these concepts is given at Fig. [0 

The usage of the directed incidence matrix allows us to rewrite model ([T]) in a vector 
form as follows: 

6 = — i diag(ljv — eZ 2 ^J B sin (^B T 6'j , (4) 

where Z 2 := (Zf ,..., Zjf), ljv := (1,..., 1) G and diag(.) stands for the matrix with the 
elements of a vector on the leading diagonal, and 0 elsewhere. 

The square of the global order parameter can be expressed as 

« 2 = + 

\ 3<k 

However, to build our Lyapunov function, we define a generalized norm of order r as 



r 


2 


E- lJcos(R T 0) 
N* 


( 5 ) 
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Note that R 2 requires the sum of all cos (9j — 9k) with j < k (for j,k = 1,..., N), but 
its generalization r 2 takes the sum (lj cos(f? T d)) only through the edges of the graph. In 
the case of full coupling graph, direct substitution yields that both global and generalized 
norm of the order parameter have the same expression . 

For any connected symmetrical coupling graph, one can check that the maximum of r 2 
is the unit, and that R 2 = 1 if and only if this value is achieved 23]. 

Let 

U{9) = 1 - r 2 (6) 


be a candidate Lyapunov function. It is clear that the minimum value of U(6) = 0 cor¬ 
responds to the maximum value of r 2 = 1, which is equivalent to the fully synchronized 
state. 

In fact, algebraic manipulations reveal that 


U(9) 



and that the differential of U is given by 


DU=^(Bsin(B T 9)) T . 


( 7 ) 

( 8 ) 


As a result, we synthesize in the next theorem the previous intuitively stated idea that 
if £ is small enough then full synchronization is a robust phenomenon related to small 
perturbations over initial conditions. 


Theorem 1 . In Model (Qp, if e is smaller than a critical value e c := l /d 2 ma:r , then the syn¬ 
chronized stated (R— 1) is Lyapunov stable. 


Proof. Consider the potential field U(9) defined in ([6]). So, using the vector form of the 
model 03]) and the expression of the differential D U from OS]); we have that ^U(9{t)) equals 
to 


2 N 2 


(sin (^B t 9^ L? T diag(ljv — eZ 2 ^) B sin(f? T o') , 


If we set x := B sm(^B T 9^) , then we have that x T diag(ljv — eZ 2 ) x is larger or equal 
than (1 — £d 2 iax )||x|| 2 . Moreover, we can also define a lower bound for ||a;|| 2 , since ||x|| 2 = 
sin(L> T 6^ B t B sin(i? T O') > A 2 (-B T -E>) sin(i? T 6^ = 2X 2 (L) sin {p T O^j ; where A 2 (L) is 
the algebraic connectivity of the graph. In the last inequality we used that 1/2 BB T = L and 
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that both matrices BB T and B 1 B have the same non-trivial eigenvalues 0 < A 2 (L) < ... < 
A a r(L), where A 2 (L) is strictly larger than zero because the coupling graph A is connected 


24], Therefore, 


d U(e(t)) < -4jA 2 (L)(l-eCd|| s in(fl T 0) 


d t 


As a result, e < e c := ’/d? riax implies that f^U(9(t)) < 0, then the fully synchronized state 


R = 1 is stable. 


□ 


IV. DYNAMICS OF PARTIALLY SYNCHRONOUS STATES 

In this section numerical simulations are performed to illustrate the rich repertoire of 
behaviors that model (JT]) may exhibit, specially beyond the threshold £ > £ c , where Theorem 
|T] cannot be applied. 


A. Quantification of dynamical regimes 


The numerical integration scheme applied is a forth order Adams-Bashforth-Moulton 
Method (see 25]) with discretization time step h = 0.01. We calculate the partial synchro¬ 


nization metric s from Ref. 


26], which, for every two oscillators i,j in the network, takes 


values Sij(I) G [0,1], indicating how much the mean phase difference between and 9j varies 
in the time interval / := [ti,^], with t\ < £ 2 - This metric is defined as 




ft 2 


^2 t\ Jti 


;i (^ 


One can check that if 9ft) = 9j(t) + p for some constant p, then the exponent in the 
previous integral is constant and %(/) = 1. Nevertheless, if Oft) — 9j(t) mod 2ir assumes 
every possible value over the unit circumference with not clear trend, then Sij(I) is close to 
zero. Now, we average contributions of all neighbor oscillators i, j under a graph A with N 
nodes to write 


N 


^ ij =1 


where E is the quantity of undirected edges in the graph. 

To exclude transients and to detect the statistically stationary state, we adopted the 
following procedure. For all experiments the time interval [0, 2.10 3 ] is always considered as 
transient time. Then, the numerical integration is performed in the subsequent intervals 
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h [(k — l),/c] 10 3 , with k > 3, until the first k = k such that — s(I *.) < 0.01, 

or k = 10 is achieved. Only such a time interval Ij, is regarded as non-transient. For 
the subsequent analysis, we use values of the phases 6{t) in the stationary time interval 
regime Ij, (whose beginning is shifted to t = 0 without loss of generality) at points t E I : = 
{fh,ie{0,l,...,10 5 -l,10 5 }}. 


B. Examples of complex behaviors 

As it was claimed before, in dependence on the network structure, very different types 
of the dynamics are possible. In order to give impression on it, we present simulations of 
model dU) with two different coupling graphs displayed as inserts in Fig. [2j Both networks 
have N = 10 nodes and they differ only by the rewiring of a single edge. We performed 
simulations for 10 random initial conditions chosen with uniform distribution over [0, 2i r] 
for each experiment. For all these initial conditions l = 1,..., 10, the norm of the order 
parameter R l (t) , according to (JUT) , is computed from the time series. As explained in the 
previous section, in these calculations a transient time is eliminated and that a statisti¬ 
cally stationary regime / of 10 3 units of time and := 10 5 + 1 points is considered. 
Then, also for each distinct initial condition, the maximum, average and minimum val¬ 
ues of the associated norm of the order parameter are computed, respectively denoted by 
ALax : = ma x teI R l (t)] (R l ) := (#/) 1 'Etei Rl ( t )> and R Lm ■= min tG jR l (t). Of course, 
R l (t) converges to a constant if and only if R l max = (R l ) = R l min . Now, having different 
simulations for a fixed coupling graph, we evaluated the maximum, average and minimum 
value of the average value of the norm of the order parameters over this ensemble, respec¬ 
tively denoted by max{(A)} := max; = i r ..,io(-R < ); mean{(A)} := (10) 1 T,i=i,...,io(R l ); and 
min{(i?)} := miip =1) 10 (A i ). So, if the norm of the order parameter converges to the same 
value for all initial conditions simulated, then max{(i?)} = mean {(A) } = min{(A)}. For 
the cases where the norm of the order parameter does not converge over all initial con¬ 
ditions, it will be useful to examine the overall maximum and overall minimum values of 

the norm of the order parameter, respectively denoted by max{A max } := max; = i.i 0 R l max ] 

and min {A min } := miip = i.io R-mm- Thus, if there is no fixed phase synchronization for all 

the initial conditions simulated, but the norm of the order parameter presents only small 
deviations around a common value, then the gap between max{A max } and min{A m i n } is also 
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FIG. 2. Numerical results for Model (HD as a function of e, for the coupling graphs despited as 
insect, including 10 random initial conditions. A black line corresponds to mean{(i?)}, while the 
interval between min {(/?)} and max{(i?)} is shown as a gray strip. The gap between min{i2 m ; n } 
and max{i? max } is shown as an orange strip. Since the orange strip is by construction larger or 
equal than the gray one, the first one is not displayed in the figure when they coincide. Left 
vertical axes show values related to norm of the order parameter, while the right ones represents 
the maximum Lyapunov exponent A max , shown as a red dashed line. Letters in green vertical lines 
from the upper experiment correspond to subfigures in Fig. [3j 


small. Also notice that miri{/?, min } < min{(i?)} < mean{(i?)} < nnix{(/i')} < max{fl „ i ;it }, 
since fq nin < (R 1 ) < R l max for all initial conditions. Finally, the maximum Lyapunov expo¬ 
nent X l max for each initial condition is also computed, according to the algorithm in 27]. The 
maximum Lyapunov exponent over all the chosen initial conditions X l max max; =li 10 
is also analyzed. 


We now describe different regimes observed in the networks, using also Fig. [31 where we 
depict time series of R(9(t)) for some particular choices of e, indicated as green letters in the 
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FIG. 3. Evolution of R(t ) for different values of e indicated in green at the upper experiment from 
Fig. [2 Every color represents a different initial condition, while pairs of solid/dashed lines with 
the same color correspond to solutions whose initial conditions differ not more than 10 -4 at each 
coordinate, (a) e = 0.04: full synchronization; (b) e = 0.08: fixed phase synchronization; (c,d,e) 
e = 0.15,0.28,0.35 respect.: examples of multi stability; (f) e = 0.70: example with A max > 0. 

upper panel from Fig. [2] (this is the case we choose for illustrating different regimes). Notice 
that d max = 4 in both cases, so Theorem [T| guarantees that for e < e c = 1/4 2 = 0.0625 
the full synchronization state, R —» 1, is locally stable as illustrated in Fig. [3] (a) (with 
£ = 0.04). 

Panel (a) in Fig. [3] illustrates full synchronization in the network for e < e c . For e slightly 
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FIG. 4. Example of group formation. Details of one of the trajectories from Fig. [3] (b) e = 0.08. 
On the left side, the coupling graph with s(i,j ) in its edges is shown. On the right side, a histogram 
of 6{ — ip is presented with color code representing the normalized frequency. 



bigger than £ c , simulations suggest that a stationary regime of partial phase synchronization, 
where R —> c < 1, is locally stable as shown in Fig. [3](b) (e = 0.08). Details of this state 
are clear from Fig. El There we show the that the synchronization between the individual 
oscillators is complete if measured by quantity ,s t j , and all the oscillators have the same 
frequency. However, the oscillators are split into two groups with a constant phase shift 
between them; this division originates in the edge which connects the two largest hubs in 
the network (vertexes 1,8). 

For larger values of £, the regimes are still static but with multistability. For instance, 
at e = 0.15 (see Fig. [ 3 ] (c)) two stable configurations emerge with R —> c, with c ~ 0.471 
(black) ore ~ 0.511 (blue), depending on the initial condition. Fig. [5], which is analogous to 
Fig. El shows the existence of three subgroups, whose members may vary according to the 
initial condition. 

Other types of multistabilities appear for instance at £ = 0.28 and £ = 0.35, as illustrate 
in Figs. El (d,e). For £ = 0.28 (panel d) some initial conditions do no converge to a fixed 
phase synchronization, but to a regime where the order parameter R is periodic in time. For 
£ = 0.35 (panel c), the norm of the order parameter of all trajectories simulated becomes 
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FIG. 5. Example multi-stability with group formation. Details of two trajectories from Fig. [3] (c) 
e = 0.15, the left picture corresponds to the solid black line and the right one to the solid blue. 
Histograms of 6i — ip are presented with color code representing the normalized frequency. 

periodic. Fig. [5] provides an illustration of this regime. 

Finally, for e = 0.70 (Fig. [3] (f)), one observes a chaotic state with A max > 0, the 
distribution of phases and frequencies is illustrated in Fig. [7] 

If £ G [1,1.5], we also obtained multistability, with the coexistence of solutions converging 
to phase-lock and irregular order parameter after the transient, similar to Fig. 0[d). 

Now, we compare the results for two slightly different networks depicted in panels (a) 
and (b) in Fig. 3 The interval of values of e with fixed phase synchronization for all 
initial conditions simulated is very similar for both networks, namely e c < e < 0.25; also 
multistability of static partial synchronous regimes have been observed in both cases. 

When e G [1,1.5], contrary to case (a), we observed that the solution for all initial 
conditions converged to the same phase-lock regime, similar to Fig. [3] (b). 

In conclusion of this section, in Fig. Owe show simulation results for two other networks. 
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sin(0i-i/f) 


FIG. 6. Example of periodic norm of the order parameter. Details of one of the trajectories from 
Fig. [3] (d) e = 0.28. On the left side, the coupling graph with s(i,j) in its edges is shown. On the 
right side, a histogram of 0* — ^ is presented with color code representing the normalized frequency. 
We denote by ip{t ) the argument of the order parameter. The bottom picture shows that the curve 
(sin(#i(f) — ,sin(0g(f) — ip(t))) is closed. 


Panel (a) shows a random network with N — 10 nodes and 20 undirected edges. Here 
predominantly static regimes are observed, only in small ranges of coupling constant chaos 
with a positive Lyapunov exponent appears. Static regimes, however, demonstrate a large 
degree of multistability. In panel (b) we show a scale-free network with N = 50 nodes and 
100 undirected edges. Here static states are rare, typically irregular regimes with low values 
of the order parameter are observed. 


C. Dependence of partial synchronization regimes on network structure 

We have seen that partially synchronous states can be rather different even for similar 
networks. It is therefore difficult to make general predictions for a relation between the 
network properties and the dynamical behaviors. Here we attempt such a description, 
focusing on the property of abundance of static regimes in comparison to time-dependent 
ones. For this purpose we define the convergence index I c as the ratio of values of e e [0,1.5] 
such that R converges to a constant value, considering all the 10 random initial conditions. 
So, both networks Fig. [2] have similar values of the index: I c « 0.530 in case (a) while 
I c ~ 0.549 in case (b). In contradistinction, network shown in Fig. [S](a) has very large value 
of the index I c & 0.946, while that in Fig. [S](b) a rather low value I c « 0.064. 

In order to explore which features of the coupling graph are related with / c , we performed 
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FIG. 7. Example of trajectory with A max > 0. Details of one of the trajectories from Fig. [3] (f) 
e = 0.70. On the left side, the coupling graph with s(i,j ) in its edges is shown. On the right side, 
a histogram of (9,; — -0 is presented with color code representing the normalized frequency. 


numerical experiments with three sets of graphs, with N = 10, 50,100 nodes. Each set 
consists in three common types of networks, each one with 10 members, generated as: (i) 
random (Erdos-Renyi) graphs with 2N edges; (ii) scale-free graphs, also with 2IV; and (iii) 
tree graphs (N edges). The Barabasi-Albert algorithm is applied for the last two types of 
networks (ii),(iii), with an initial clique of mo nodes and with other nodes been connected 
to m existing ones. For the 21V-edges scale-free graphs, we fixed m o = 5 and m = 2; while 
for the tree graphs (N edges scale-free graphs), m 0 — m — 1. We point out that all graphs 
created are connected and symmetrical. Additionally, three sets of 10 initial conditions 
do G R^, with uniform distribution over [0, 27 t] and N = 10, 50,100, have been explored. So, 
for each of the 90 coupling graphs we computed its correspondent I c values by numerical 
integration of model P) for e = 0, 0.01,..., 1.49,1.50. 

In Table [□ we report the mean value and the standard deviation of I c for each topology 
and size of coupling graph. From these data we see that the mean value of I c increases if 
we go from tree to scale-free and to random graphs, respectively. However, these difference 
becomes less noticeable for larger values of N. Both the mean value and the standard 
deviation of I c decrease with larger networks. 
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FIG. 8. Numerical results for Model © as a function of e, for the coupling graphs depicted as 
insect, including 10 random initial conditions. The legend of the pictures is the same as in Fig. [2] 


Network 

N = 10 

II 

Ol 

O 

N = 100 

Tree 

0.421 (0.260) 

0.016 (0.006) 

0.008 (0.004) 

Scale-free 

0.857 (0.029) 

0.050 (0.015) 

0.013 (0.003) 

Random 

0.872 (0.090) 

0.183 (0.063) 

0.077 (0.022) 


TABLE I. Mean value of I c and its standard deviation (in brackets) for each network type and size 
simulated. 


We have explored different networks metrics, searching for one mostly correlated with 
the convergence index I c . Let 0 = 71 < 72 < ... 7,/v denote the Laplacian eigenvalues of 
the coupling graph 28]. Recall that this graph is assumed to be simple and connected. We 


stress that these eigenvalues express fundamental characteristics of the graph. For instance, 
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FIG. 9. Convergence index I c versus 7 *. Networks with N = 10,50,100 nodes are represented as 
circles, squares and triangles, respectively. The two experiments from Fig. [2] are shown as disks. 
We show in red an exponential fit f(x) = e 1 -' 6 ' 6 ~ L0894x for the data. 

72 is related with graph diameter and 77 V with its largest degree size. 

We found that the quantity 7 *, defined as the ratio between the maximum eigenvalue 
and the average of the non-trivial eigenvalues of the Laplacian matrix of the graph, is rather 
suitable for this purpose. Formally, it is defined as 

fori*)"' 

In Fig. |a correlation plot between I c and measure 7 * for the correspondent graph is 
presented. From there, we observe a clear trend indicating that the greater the value of 7 * is, 
the smaller is the value of I c . Independently of the network type and size, static regimes of 
partial synchronization, full synchronization and phase-lock, are typical for values of 7 * <3, 
like in the experiments from Fig. [21 On the other hand, graphs with larger values of this 
measure yields more irregular dynamics, like time-dependent periodic and chaotic regimes, 
as the ones from Fig. [HI 

V. CONCLUSION 

In this work we introduced and studied a Kuramoto-like model of identical oscillators with 
non-linear coupling. Our main parameter was e, which governs the coupling nonlinearity 
strength. It is clear that the most influence of nonlinearity in the coupling is on the hubs 
which experience strong forcing from many connected oscillators, while less connected nodes 
may still operate in a linear-coupling regime. 
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We proved that if this parameter is smaller than the inverse of the square of the maxi¬ 
mum vertex degree in the network, then the full synchronized state is stable. Via numerical 
experiments, we showed that our model can display a variety of other qualitative behaviors 
of partial synchronization, like stationary phase locking, multistability, periodic order pa¬ 
rameter variations, and chaotic regimes. We explored the relative abundance of stationary 
phase locking regimes under different network topologies. Statistical analysis performed sug¬ 
gests that tree graphs are much less likely to exhibit stationary phase locking in comparison 
with scale-free or random networks. In addition, this type of behavior becomes more rare 
if we increase network sizes, irrespective to the network topology. Finally, we also found a 
good correlation between the ration between the maximum eigenvalue and the average of 
the non-trivial eigenvalues of the Laplacian matrix of the graph, and the proportion of the 
repulsion parameter values which yield stationary phase locking. The greater this measure 
is, the smaller tend to be presence of stationary phase locking states in the system. 

As a future research, we want to investigate analytical conditions and correlations involv¬ 
ing other graph measures related to other forms of synchronization in the model. 
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